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We study the entanglement entropy and spectra of a coupled array of N one dimensional quantum 
Ising chains in their continuum limit. Employing a DMRG algorithm specifically adapted to the 
study of coupled, continuum systems, we are able to study large arrays of chains (up to N=200) both 
in their gapped phase and in the approach to criticality. Away from criticality the entanglement 
entropy obeys an area law. Close to criticality the entanglement entropy continues to obey the 
area law but possesses an additive piece scaling as c e // log(A r )/6 with c e // ~ 1. We also study the 
entanglement spectra of the coupled chains. Away from criticality in the disordered phase the low 
lying portion of the entanglement spectra appears similar to that of a single gapped quantum Ising 
chain. As the critical point is approached the entanglement gap closes. A finite scaling analysis 
shows that the entanglement gap and the energy gap vanish at the same value of interchain coupling. 

PACS numbers: 



In the past decade concepts borrowed from informa- 
tion theory have become important tools in analyzing 
the properties of low dimensional strongly correlated con- 
densed matter systems pQ. A particularly useful quan- 
tity in this regard is the bipartite entanglement entropy, 
Se- Se, a measure of non-local quantum entanglement 
in a system can be used to characterize quantum critical 
points PHI], access whether a given model has hidden 
topological order [5HB], and provides a simple measure 
by which to judge whether a density matrix renormal- 
ization group (DMRG) algorithm, one of the most com- 
monplace numerical techniques in low dimensions, will 
succeed [3j|9]. 

Most is known about Se in one spatial dimension (ID). 
In ID Se signals the onset of criticality through an as- 
sociated universal logarithmic divergence in the system 
size, L plH]. Both the coefficient of this divergence (i.e. 
the central charge of the theory's conformal algebra) to- 
gether with its subleading corrections in L (determining 
the theory's operator content [101 E]) serve to uniquely 
specify the underlying critical theory. 

In spatial dimensions greater than one, less is known 
on how the behavior of Se reflects criticality. Regardless 
of criticality, Se, possesses a term scaling as the area of 
the boundary separating the bipartite region 12J. Be- 
yond the area law term, there can be subleading, univer- 
sal contributions to Se- Generalizing the results in one 
dimension, the AdS/CFT correspondence suggests that 
Se in all odd spatial dimensions will be characterized by 
universal logs [T31 [H] . Universal terms have also been ar- 
gued to appear in a set of critical theories in 2D known as 
conformal quantum critical points (CQCP) [3 IT5H20] as 
well as systems with spontaneously broken symmetries 
[2"Tl |2"2"] . Recent studies of gapless states on the torus 
[2"0l |2"51 12~3] have confirmed the existence of apparently 
universal terms that depend on system size and shape. 

In this letter we investigate the presence of universal 
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FIG. 1: a) An array of continuum chains of length, R, coupled 
together with an interchain coupling, J±. Each chain is a 
single site in our DMRG algorithm, b) The phase diagram of 
coupled QICs. 

terms in Se for 2D critical systems. We do so by studying 
these systems in a particular limit, one where the model 
is represented as a mixture of continuum and discrete 
degrees of freedom. We consider as a particular exam- 
ple the two dimensional quantum Ising model. We study 
such a mixed system because it is amenable to treatment 
by a ID-like DMRG algorithm, as will be explained. We 
find that the 2D entanglement entropy shares a number 
of characteristics with systems in ID, in particular, uni- 
versal logarithmic scaling. 

With our DMRG algorithm we are not only able to 
study Se but also the entanglement spectra (ES). While 
the entanglement spectra was initially introduced as a 
way of detecting topological order PS] , it has been stud- 
ied in non-topological systems in order to attempt to elu- 
cidate connections between ES and the system's ordinary 
excitation spectra [2BH22|. Here we will show that the 
entanglement spectra can be used a robust means to de- 
tect a 2D critical point inasmuch as the entanglement gap 
vanishes at the critical point. We will show that the scal- 
ing of this gap, like Se, bas ID-like features as described 
in Ref. [301 E]- 

Model and DMRG algorithm: We study the 2D 
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quantum Ising model where we envision the model as 
a set of ID quantum Ising chains (QICs) of length R 
coupled through their spin operator: 

H = Y, H l DQI + J ±Y, [ dxa^ajix), (1) 

where i is a sum over chain index. The Hamiltonian 
Hj D ® l is taken in its continuum limit, that of a massive 
Majorana fermion, H\ D ^ 7 = J dx{i^)d x ^} — iipd x t[)) + 
iAipip) where ip/tj> are right and left moving components 
of the Majorana fermion. In lattice notation, H 1D — 

—JJ2j a j a j+i + (.9 + an d we identify a z — > a, A = 
gj. In this letter we will focus on systems built from 
chains with negative mass (A < 0), as it is in this case 
that one can approach criticality - driving the system 
through a phase transition by increasing the magnitude 
of the interchain coupling, J± (see Fig. [I]). To study 
this model we employ a density matrix renormalization 
group (DMRG) algorithm adapted to coupled ID chains 
as described in [32]. Like all DMRG algorithm's this 
allows us to extract readily the block Se and spectra. In 
this realization of DMRG, we treat individual chains as 
equivalent to individual lattice sites in the conventional 
DMRG algorithm. 

This methodology is based in part on the truncated 
spectrum approach (TSA) to studying perturbed confor- 
mal and integrable field theories [33] ■ I n this method 
the underlying conformal or integrable theory provides a 
particularly apt basis in which to study relevant (in the 
renormalization group sense) perturbations. With a rele- 
vant perturbation, the low energy sector of the full theory 
can be understood as a mixing of the low energy sector 
of the unperturbed theory (even if its energy spectrum 
is dramatically different). Thus the high energy part of 
the theory need not be considered or can be taken into 
account in a variational scheme borrowed from the nu- 
merical renormalization group (NRG) [34) . Thus in im- 
plementing this chain DMRG we place a cutoff in energy, 
E c , on the uncoupled chain Hilbert space. 

Like perturbations in the TSA, this DMRG trades on 
being able to compute matrix elements of the interchain 
coupling exactly, i.e. (s\<jj<jj+i\s'), where \s), \s'} are two 
states on a pair of neighbouring chains. The ability to 
compute these matrix elements means we are able to in- 
corporate much of the strongly correlated physics before 
the numerical analysis has even begun. 

The standard DMRG algorithms in two dimensions 
run into difficulties because Se grows with the length of 
boundary between blocks (in this implementation, R, the 
chain length) . Here our use of continuum chains plays an 
important role in allowing the DMRG algorithm to work 
successfully. In continuum field theories, the finite size 
errors are exponentially suppressed in system size, R [35]. 
This means that we can keep R comparatively small but 
still be in the effective infinite volume limit. And if R 
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FIG. 2: Se as a function of the number of chains, N, in the 
disordered phase and near criticality. Inset: Se for the same 
phase as a function of R, the chain length. 



is small, the corresponding block Se can also be kept 
small allowing the algorithm to operate successfully. As 
a corollary to this we find we need to keep comparatively 
few eigenstates of the reduced density matrix, ranging 
from the 10's deep in the ordered phase to at most 120 
close to the critical point in order to obtain truncation er- 
rors on the order of 10 -5 (for additional details see 37 ). 
Thus while the Hilbert space of the individual sites (i.e. 
the chains) can involve many hundreds of states, we need 
to keep far fewer states from the reduced density matrix. 

Because of these features, this DMRG algorithm has 
been able [35] to successfully analyze various properties of 
large arrays of coupled QICs. It was able to reproduce the 
scaling form (given in terms of the dimensionless combi- 
nation Jj^A -1 ) of the first excited gap in the disordered 
phase. For this system it is possible to analytically com- 
pute the finite chain length, R, corrections. In [32) an 
excellent match between these analytical computations 
and the DMRG numerics was found. But most signif- 
icantly, the exponent v governing how the gap, A exc , 
vanishes as the critical coupling, J± = J c , is approached, 
i.e. A exc ~ |Jj_ — J c \ u , was computed. There it was 
found v = 0.622 ± 0.019, in reasonable agreement with 
the accepted value, v = 0.630, for the 3D classical Ising 
model. This last computation shows the DMRG suc- 
cessfully captures 2D physics despite using ID chains as 
building blocks. 

Entanglement Entropy: Starting from disordered 
chains, A < 0, and coupling them together with in- 
creasing Jl, we study the behavior of Se and the en- 
tanglement spectra, both deep in the disordered phase 
and in the vicinity of the critical point. In Fig. [2] we plot 
some typical results for Se in the disordered phase with 
R = 10. At very small J±, Se is essentially constant with 
N - saturating even for systems of fewer than 10 chains 
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FIG. 3: Se close to criticality for a block of length N/2. 
Left: as a function of the chain circumference, R. Right: as 
a function of total length, N. 



- but increases linearly with R, Se ~ R, hence display- 
ing an area law in accordance with expectations. On 
increasing J± the value of N at which Se saturates also 
increases and for J± > 0.18 a logarithmic dependence 
on N becomes evident. At J± = 0.195 this Se ~ log A" 
behavior extends up to systems of 7Y ~ 100 chains. Note 
that at this coupling Se still obeys an area law (Fig. |3| . 
For J± > 0.196 the value of N at which Se saturates 
begins to decrease again. Our findings are compatible 
with a scenario in which the system passes through a 
critical point while Se tracks the correlation length £: 
when J± <C J c the system is deep in its disordered phase 
and £ -C Na (here a is the interchain spacing which we 
set equal to unity) with only local entanglement so that 
Se is independent of N. As J± approaches J c , increas- 
ingly long ranged correlations build non-local entangle- 
ment into the system and Se scales as some non-trivial 
function of N. Once J± exceeds J c the system enters the 
ordered phase and correlations diminish again, leaving 
Se increasingly independent of N. Because of finite size 
effects we expect the true value of J c to differ somewhat 
from 0.195. We return to this issue later. 

In Ref. [24] a Quantum Monte Carlo (QMC) study re- 
vealed intriguing ID-like 'chord scaling' of Se for several 
gapless 2D quantum ground states on the torus (in par- 
ticular the spin- 1/2 Neel state). By varying the partition 
(block) length x on a torus of total length N the authors 
of that work found that Se depended on x via a term 



C(N) log [sin (— J 



(2) 



A somewhat similar log term was also recently found for 
a CQCP in Ref. [20]. In our DMRG calculations x cor- 
responds to the number of chains in the 'system' block 
while there are N — x chains in the 'environment' block 
(Fig. [IJ. In Fig. [4] we plot Se versus the expression 



0.9 



o.s 



se R=10, N=20 




mR=10, N=40 




<m>R=10, N=100 




XXR=12 N=20 




- ^R=12'n=40 




- A=-l, J ± =0.195 


^^^^^^^^^^^^^^^^^^^ 







^2 -1.5 

ln(sin(Jtx/N ) 



FIG. 4: Se for J± = 0.195 as a function of block length. 



in Eq. [2] at J± = 0.195 for a variety of aspect ratios; 
there is a clear linear relation (though apparently the 
system is larger than the correlation length for N = 100 
at Jl = 0.195). In Ref. [23] a strong dependence of the 
gradient C on N was found, though it appeared to be ap- 
proaching a constant value for larger N. Our calculations 
for larger system sizes display only a weak dependence on 
N with C = 0.140(5) (see Table [TJ we only fit the linear 
piece of the N — 100 data). 

In ID for critical models with open boundary condi- 
tions Se takes the form 
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where c is the conformal charge. If we draw a com- 
parison with this case then we expect that at critical- 
ity the entropy of our 2D system should also scale as 
log/Y when x = N/2, with the same coefficient, C. 

We find S E = C log N -\ , at fixed x = N/2 with 

C = 0.155(5) (see the right panel of Fig. ||. This close 
agreement between C and C is suggestive of the defini- 
tion C = C = c e ff/6, with c e ff w 1. Note that our 
system is not strictly critical: we work at finite size and 
there is still a small gap at J± = 0.195, so we expect some 
departure in C and C from the true thermodynamic val- 
ues. 

TABLE I: Values of C from data in Fig. [1] 
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Entanglement Spectra: In principle much more in- 
formation is encoded in the full spectrum of the re- 
duced density matrix, p r , than in the number Se = 
—Trp r log p r alone. The ES is usually defined in terms of 
a fictitious Hamiltonian, p r = exp(— He) so that the ES 
'energies' are w = — logp r . We first consider the entan- 
glement 'gap', A^s, the difference between the two low- 
est lying values of the ES. Calabrese and Lefevre have 
shown that in a ID conformal system of finite length 
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FIG. 5: Finite size scaling of both (left) the entanglement 
gap and (right) the real gap for DMRG calculations in which 
eigenvalues of p r are kept if they are > 5 x 10 -7 . 



N, A ES ~ const/ log(7V/7r) [SU] EJ- In Fig. ^ we plot 
Asslog(AT/7r) against J± for a variety of systems with 
the same aspect ratio (R/N) but different sizes. Tellingly, 
the curves cross at a single point J c = 0.186(2) indicat- 
ing that we can use this finite size scaling of the entangle- 
ment spectrum to discern the critical point of the system. 
To support this claim we also perform finite size scaling 
on the true energy gap (E\ — Eq) in our system (right 
panel of Fig. [5| and find J c — 0.185(2). This also agrees 
well with the RG improved value J c = 0.184(3) found 
in Ref. [33] Calculating E\ — Eq is considerably more 
difficult than Aes as the former requires targeting the 
first excited state with the DMRG algorithm. This state 
has larger Se than the ground state and hence carries a 
greater numerical burden. 

Finally we consider the ES as a function of the momen- 
tum, fc, along the chain direction. It has been shown that 
the ES of spin ladders closely resembles the true energy 
spectrum of a single spin chain [2"6"H2"5] , The spectrum 
of the QIC separates into two sectors, Neveu-Schwarz 
(NS) and Ramond (R) [37]. For a disordered chain these 
correspond to even and odd numbers of solitons along 
the chain respectively. Similarly the ES splits into two 
sectors, depending on whether the state has an even or 
odd number of chains in the NS sector (assuming N/2 is 
even). Fig. [6] shows that at J± = 0.13, far from critical- 
ity and where short range entanglement at the boundary 
dominates, the low lying ES resembles that of a single 
QIC where the one and two soliton sectors are mimicked 
by ES states with odd and even numbers of NS chains. 
Closer to criticality, at J = 0.18, the ES does not re- 
semble that of a disordered QIC, in particular Aes — > 0. 
A perturbative calculation for weak intrachain coupling 
gives uj = 21og(A 2 + k 2 ) + const for the lowest 'band' in 
the ES (see [37]). There is good agreement between this 
prediction (with A = —1) and the Jj_ = 0.13 spectrum 
as shown in Fig. [6] 
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FIG. 6: Plot of the entanglement spectra both away from 
criticality (left), J± — 0.13 and closer to criticality (right), 
J± = 0.18. iV = 160. In both cases the ES is measured 
relative to its lowest eigenvalue. The two sectors of the energy 
spectrum of a QIC are also plotted, rescaled so that the lowest 
band overlaps the lowest ES band at k — 0. The curve labelled 
'theory' is the perturbative calculation in [37] . 



The ID like features that we see in our 2D system sug- 
gests the following interpretation. Using the intuition 
that comes from our treatment of a 2D system as a mix- 
ture of discrete and continuum degrees of freedom, any 
2D system can be thought of as a set of coupled con- 
tinuum chains. At a critical point, the representation of 
the system as such a mixture should not affect the criti- 
cal properties (provided the critical point is a point and 
not a line where a lattice vs continuum treatment might 
control where along the line one ends up). If at the 2D 
critical point, a finite number of chains become critical 
with the remaining chains massive with a gap of at least 
A m i n , one would expect to see ID scaling. 

In summary, we have shown that an unconventional 
DMRG technique can be used to study the entanglement 
content of 2D quantum systems. Using this technique we 
have established the existence of an additive logarithmic 
piece in Se with a possibly universal coefficient c e //, ~ 1 
for the 2D quantum Ising model. We have also shown 
that the ES gap can be used as an efficient tool to find 
the critical point. 
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Supplementary Material 

Details of Numerics: We consider arrays of N chains 
with N < 200 and A = — 1 (which sets the overall energy 
scale). The primary limitation on our numerics is storage 
for the states of two chains, necessary for the first basic 
step of the DMRG where we consider a system of four 
chains. The number of these states grows rapidly with 
energy cutoff and chain length, R. In the following we set 
E c = 8.0 and R < 12 corresponding to the incorporation 
of < 18000 two chain states. 

On approaching the critical point Se increases. This 
indicates that the spectrum of the reduced density ma- 
trix is no longer dominated by just a small number of 
eigenvalues. Consequently the DMRG algorithm must 
retain a larger number of eigenvalues in order to accu- 
rately represent the target state with a larger computa- 
tional penalty. In our DMRG algorithm whether or not 
an eigenvalue is kept in each step is controlled by the con- 
dition p re d.i > p c where p c is a threshold value. For the 
detailed plots of the energy gap, E\— E , and Aes close 
to criticality (Fig. [5]) we found it necessary to decrease 
p c by a factor of 5 (effectively increasing the number of 
kept states by 50%) to produce sufficiently accurate re- 
sults. In Fig. [7] we plot data obtained using the larger 
(less accurate) value of p c . The effect is most clearly 
visible for the larger system sizes and especially for the 
plot of the energy gap, because the first excited state 
carries more entanglement relative to the ground state 
and is thus more challenging for the DMRG to represent, 
reducing the accuracy and rate of convergence. It was 
necessary to increase the number of convergence sweeps 
in the finite volume part of the DMRG routine from ~ 4 
for the ground state with R — 6, J± = 0.17 to ~ 20 for 
the first excited with R = 9.2, J± = 0.185. We also found 
it necessary to choose our system sizes such that N/2 was 
even, to avoid small even-odd effects associated with the 
number of chains in a block. In Fig. [7] it is much harder 
to discern the crossing points, leading to crude estimates 
of J c as 0.187(3) and 0.185(5) from finite size scaling of 
Aes and E\ — E respectively. We do see however from 
this computation that as we increase the accuracy of the 
DMRG computation (as measured by kept states), the 
crossing points from the finite size scaling of the ES and 
the gap to the first excited state move closer together. 

Entanglement Spectra of Weakly Coupled 
Chains: Consider a system of N disordered (A < 0) 
QICs, where N/2 is an even number, with weak in- 
trachain exchange J±. The Hilbert space of a single 
QIC on a periodic interval of finite length R (as 
summarized succinctly in |35j ) splits into two sectors, 
known as Neveu-Schwarz (NS) and Ramond (R). For a 
disordered chain the states in the NS sector consist of 
even numbers of solitons with their momenta quantized 
as ki = 2'Kmi/R with mj taking half-integer values. 
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FIG. 7: Finite size scaling of (left) the real gap and (right) 
the entanglement gap for DMRG calculations in which eigen- 
values of p r are kept if they are > 25 x 10~ 7 . 



In contrast the states of the R sector consist of odd 
numbers of solitons with ki = 27m i /i? where rij takes 
integer values. The energy of a p soliton state is the sum 
J2i=i \/A 2 + k\ with p even and odd for the NS and R 
sectors respectively (we neglect the two vacuum state 
energies, the difference between which is exponentially 
small in the length R). 

The ground state of the N chain system for J± = is 
then a tensor product of zero soliton (vacuum) NS states. 
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i^o) = n i ns >, 



(4) 



Treating the intrachain coupling in the Hamiltonian as a 
perturbation 



H± = J± dr S~] <7*(r)er <+ i(r), 
J o i=l 



(5) 



and using translational invariance along the chains, we 
find the first order contribution to the ground state is 



IV>i 



N-l 

=-*.*e( n i ns ), 
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|fc)J-fc) m (fcHNS) (-fcHNS) 
2%/A 2 + k 2 



(6) 



where we have used the property that matrix elements 
of the spin operator only connect different sectors. Fur- 
thermore we have assumed that the most important con- 
tribution comes from the one soliton (Ramond) states, 
\k), as there is a mass gap, A. For large R the matrix 
elements simplify [35] and we obtain 



J -^"t( II in%)E" : 
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2A 



+ 1 



A 2 + k 2 



(7) 
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where a = x 1.35783834 ■ 

corrected state is 

N = (V'ol^o) + 



The norm of this 



J ± a 2 



L2A(A 2 + fc 2 ) 



(8) 



As the normalization is 0( Jj), the same order as the low 
lying entanglement spectra, it is necessary to take this 
normalization into account and renormalize the state as 
■A/" -1 (IV'o) + l^i))- We now form the density matrix 
p = W^d'f/'o) + |^i))((^ol + (V'll) and ta ke the trace 
over the states of the N/2 chains on the right-hand side 
of the system. The resulting contributions to the reduced 
density matrix, p r , can be separated into even and odd 
sectors depending on whether they include an even or 
odd number of NS chains (equivalently an even or odd 
number of total solitons) on the leftmost N/2 chains. The 
zero soliton contribution is given by 
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(9) 



where the N dependence cancels on expanding JV^ 1 and 
keeping terms in p° to 0(J±) only. We extract the 
'ground state energy' of the ES as 
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which is of order 10~ 2 for the range of parameters we 
use. The one soliton contributions can be organised by 
chain momentum, k: 



N/2-1 



i{k)= J] INS^NSI, 



J ± a 2 



2A(A 2 + k 2 ) 



so that the lowest ES band is given by 

Wl (fc) = -21og(^)+21og(A 2 
This is the curve plotted in Fig. [6] 
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